NPS ARCHIVE 
1966 

HALLAS, H. 



THE HYBRID SIMULATION OF A 
TACTICAL WEAPON DISCRETE FILTER-CONTROLLER 

HENRY GEORGE BQUCHIER HALLAS. 





V 







Library 

U. S. Naval Poatgraduar* ->c 
Monterey, California 



DUDLEY KNOX LIBRARY 
NAVAL POSTGRADUATE SCHOOL 
MONTEREY CA 93943-5101 



This document has been approved for public 
release and sale; its distribution is unlimited# 






■ V? 



DUDLEY KNOX LIBRARY 
NAVAL POSTGRADUATE SCHOOL 
MONTEREY CA 93943-5101 



THE HYBRID SIMULATION OF A TACTICAL WEAPON 
DISCRETE FILTER-CONTROLLER 



by 



Henry George Bouchier Hallas 
Lieutenant, Royal Canadian Navy 
B.A.Sc., University of Toronto, 1959 



Submitted in partial fulfillment 
for the degree of 

MASTER OF SCIENCE IN ELECTRICAL ENGINEERING 

from the 

UNITED STATES NAVAL POSTGRADUATE SCHOOL 

May 1966 



ABSTRACT 



This hybrid simulation demonstrates the feasibility of using 
a digital computer for the realization of a discrete filter-con- 
troller to drive a 3"/50 gun mount in response to step, ramp and 
stabilization orders. The gun position designation signals are 
subject to random environmental noise and the observation of 
position is contaminated by random measurement noise. These 
noise sources are assilhned to ej&ibit a normal distribution with 
zero mean. A filter is used to predict plant position and velocity 
from a noisy observable position, and a controller is used to 
provide the desired feedback of these states for a ripple-free 
response. The simulation results indicate that there is no de- 
tectable difference in response between constant filter gains and 
adaptive filter gains for this time-invariant plant. The results 
also demonstrate that the application of good programming tech- 
niques is paramount in minimizing the timing and scaling 
limitations imposed by hybridization. 
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1 . Introduction 



In the past few years a great deal of progress has been made 
in the theoretical application of state space techniques for the 
design of sampled-data control systems. The objectives of this 
paper are to apply such theory to the design of an optimal con- 
troller and Kalman filter for the control of the 3"/50 RF Twin Gun 
Mount and to construct a hybrid simulation of the system so de- 
signed, simulating the plant on the Pace TR20 analog computer 
and representing the filter- controller by a program in the CDC 
160 digital computer. 

The analytical methods outlined in References [1], ft 
j3^J and are applied to the design of a discrete filter- 
controller with constant filter and controller gains to replace 
the conventional continuous feedback control scheme. These 
references assume that the system plant may be expressed as 
a set of linear, constant coefficient, differential equations. 
References DJ and pj outline digital computer programs that 
compute the gain matrix elements for the filter, and the feed- 
back coefficients of the controller, respectively. Reference £4j 
presents a very thorough summary of the theory as applied to 
the optimum filter-controller problem and introduces the treat- 
ment of deterministic inputs as shown in the general simulation 
schematic of Figure 1-1 . 

System operation with constant filter gains is compared to 
that with adaptive filter gains , where the latter are calculated 
on-line in the digital computer. 

2 . Description of 3"/50 Gun Mount 

The conventional gun control system uses a typical synchro- 
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amplidyne drive system as shown in Figure 2-1. Position error 
is determined by a control transformer, the error voltage is then 
demodulated, amplified, compensated and applied to the ampli- 
dyne power amplifier. The train and elevation servos are 
essentially identical. The system has a coarse control synchro 
for slewing onto the target bearing quickly and a fine control 
synchro for tracking the target smoothly. The control perform- 
ance constraints are a maximum angular velocity of 30 degrees 
per second and a maximum angular acceleration of 75 degrees 
per second squared. 

To introduce the discrete filter- controller loop, the conven- 
tional control loop was broken at the input to the amplidyne 
power amplifier and the output shaft of the d-c drive motor as 
shown in Figure 2-1. The breaking of the system at these points 
allows either the conventional compensation loop and control or 
the discrete compensation loop and control to be used by simply 
throwing a switch . 

Reference [5] gives the following mathematical descrip- 
tion of the gun mount open-loop transfer function for train and 
elevation as determined by frequency and transient response 
techniques: 

G(s) = K(s 2 + as + b) , . 

s'(s + p 1 )(s + p 2 )(s + P3) (s^ + cs + d) ' 

where K is the forward path gain, and the three simple 
poles are introduced by the characteristics of the d-c amplifier, 
amplidyne, and the drive motor. The complex zeroes and poles 
are introduced by the mechanical structure of the gun mount 
which causes an open loop resonant condition at about 8 cycles 
per second. 
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This transfer function can be approximated by the following 
second order plant whose time constant represents the combined 
characteristics of the drive motor, and the effective friction and 



inertia at the motor shaft; 
G(s) = K 



( 2 - 2 ) 



s(s + a) 

This approximate transfer function models the real system well 

0 

enough for the purposes of this paper. 

In the real gun control system, non-linearities exist be- 
cause of saturation and because of the limits placed on velocity 
and acceleration that were designed into this control system for 
protective measures. To simplify the simulation and to be able 
to apply linear state space theory the plant is assumed to be 
linear even though this approximate transfer function of 
equation 2-2 holds only within certain bounds. 



GUN MOUNT PERFORMANCE CHARACTERISTICS 



Train 


Elevation 


Transfer 


Function G(s) = 3.0 


G(s) = 6.8 


s(s + 3.5) 


s(s + 3.5) 


Uncompensated 

Zeta _ 

Greater than one 


0.668 


Natural 


Frequency , 

1 .732 radians/sec 


2 . 62 radians/sec 



TABLE 2-1 

Assuming unity feedback the train transfer function would 
have a closed loop response with a damping ratio slightly 
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greater than one and a natural frequency of 1.732 radians/ 
second. The system is over-damped and has a very low natural 
frequency. Consequently the response would be sluggish. In 
reference , personnel at the Naval Electronics Laboratory 

outline the application of a 'dominant mode 1 type of compensator 
which introduces a more desirable pole location having a damping 
ratio of 0.7 and a natural frequency of ten radians/second. Using 
a sampling rate of 0.1 seconds and this 'dominant mode' com- 
pensator they were able to provide a digital control to the 3"/50 
gun mount and to obtain a ripple free response to step and ramp 
inputs . 

3 . Gun Position Orders 

The gun, however, must in addition respond to sinusoidal 
stabilization signals of roll and pitch as well as the determinis- 
tic inputs of a step and a ramp. Furthermore these signals are 
all contaminated by noise. 

For this paper the sea motions of roll and pitch are defined 
by the following equations: 



where 

AR is roll amplitude in degrees 
AP is pitch amplitude in degrees 
T f is the roll period in seconds, T f = 10 seconds 
Tp is the pitch period in seconds , T p “ 5 seconds 
and the magnitude of the roll amplitude is assumed to be twice 
that of the pitch for each sea state. 

In the hybrid simulation the stabilization signals are 
sampled and subsequently represented by the output of a zero 



Roll =AR*sin(w r t) 
Pitch =AP*sin(w t) 



(3-1) 

(3-2) 
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order hold. It is assumed that the random noises appearing on 
the stabilization signals are Gaussian, with zero mean and 
with standard deviations that can be determined. 

In the conventional mode of operation of the gun mount , 
stabilization signals would be applied continuously. The mag- 
nitude and direction of these signals is a function of the sea 
state. They are required to keep the gun stationary with re- 
spect to a horizontal plane. When a target is designated to the 
gun, a step input causes the gun to slew to this position at max- 
imum velocity. Upon acquiring the target, a ramp input is re- 
quired for smooth tracking of the target. 

It is the intention of this study 

1) to use a filter to estimateuthe plant state variables (gun 
position and velocity), and 

2) to use a controller to provide the right combination of 
these estimated states in forming a scalar control to obtain an 
optimum response. The implementation of these ideas, as 
depicted in Figure 1-1, will be described in the succeeding 
sections . 

4 . State Space Description of the Plant 

Consider the general second order plant of equation 2-2 to 
describe the train and elevation motions of the gun. This plant 
has been assumed to be a linear, time-invariant system and 
by the methods outlined in Appendix I can be described as a 
set of first order, linear differential equations given in the 
following matrix format: 

x = F x + D u (4-1) 
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where 



and 



F = 


0 l" 


D =' 


V 


1 


0 -a 




k 

L J 



The state transition matrix can be found as a solution of 
the unforced state equation in the following manner: 
x = F x 

s x(s) - x(0) = F x(s) 
x(s) = [sI-F ]“■*■ x(0) 
where I is the identity matrix 
but 




1/s l/s(s+a) 

0 l/(s+a) 



For the continuous case, 

*(t-t 0 ) = L' 1 [sI-F] -1 = 



1 

0 



(l-€ -a (t_t o))/a 
€ _a (t— 1 0 ) 



(4-2) 



The state transition matrix for the discrete case is 



<J> = 



1 

0 



(l-€ -aT )/a 



-aT 



(4-3) 



jT,a 



where T = t-t , and T is the sampling interval 
and a a is the time constant of the plant. 

If the control u is represented as the output of a zero 
order hold and therefore is constant over the sampling interval 
then the distribution matrix can be found as a solution to the 
forced state equation for the continuous case in the following 
manner: 
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A= £ D dt^ 



o 

A = % L" 1 



A= 



5 

t. 



K ^ ((l-c _a ^ t_t ob/a)dt 

K ^ € -a(t-t 0 ) dt 

' ^ 



l/s(s+a) 




0 


1/ (s+a) _ 




K 



dt 



(4-4) 



where K is the plant gain. 

The distribution matrix for the discrete case is 

A = 



K j ((l-€‘ aT )/a)dt 



K 



0 

T 



f-aT 



dt 



(4-5) 



T,a,K 



Then the complete discrete state equation depicted graphically 
in Figure 4-1 is given by 

x(k+l) = <tx(k) + Au(k) (4-6) 

By using the digital program 'Phidel' described in 
reference [ 2 ], the transition and distribution matrices for the 
train plant shown in Figure 4-2 for a sampling interval of one 
second were found to be 





"l.OOO 


0.277 


and 


A = 


0.475 




0.000 


0.030 






0.651 



5 . The Controller 

It is desired to determine a set of feedback controller gains 
to control the plant in accordance with a selected performance 
index. Reference [ 2 ] defines one possible performance index 
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for a discrete sampled-data system as follows: 

J(N) = Min £ |"x t (k)Qx(k) + Ru 2 (k-1) 1 (5-1) 

u(k) K=M L J 

where 

Q is a nxn positive definite, symmetric, constant matrix 
N is number of stages of the process 
R is a positive constant 
u is the control 

x is a nxl vector of plant states 
For the minimum time case and for the second order plant being 
considered in this paper the cost function reduces to 

J(N) = x^ (N) + x^ (N) (5-2) 

Reference [2] goes on to develop the recursive algorithms for 
this discrete final-value control problem as follows: 

¥ = 0, A T = 0, P Q =1 

* ■ 

A t Vi a 
^ = <f + aa t 

P = p ^ 

1 i i-1 1 

where 

i = 1 to N; N the number of stages of process for 

the minimum time case is equal to the order of the system 

I is the identity matrix 
T 

A is a lxn vector of feedback controller gains 
1 

and 

and 1 ® r are matrices involved in the intermediate steps in 
the recursive solution for the feedback gains. 
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From the Fortran 60 program 'Optcon', modified to minimize 
the final states, the following feedback controller gains were 
obtained for the train plant: 

A T = [-1.203 -0.343 ] 

A listing of program 'Optcon' is given in Appendix II. Table 
II— 1 gives the definitions for the terms for this program and in 
addition suggests values of the parameters of the cost function 
equation 4-1 for the three cases of minimization of time, min- 
imization of control effort, and minimization of both time and 
control effort. 

6 . The Filter 

In the controller discussion, it is assumed that all the 
states of the dynamic system are observable states. This con- 
dition is certainly not always true. In addition, the observable 
states may be measured only at the cost of some ambiguity due 
to measurement noise. The digital filter will provide a best 
estimate of all the states by weighting the past information with 
the present observable states. This weighting is performed with 
the knowledge of the environmental and measurement noise, 
both of which, it has been assumed, have independent Gaussian 
distributions with zero mean and a standard deviation that can 
be determined in some fashion. 

The recursive equations 6-1 , 6-2 and 6-3 for the Kalman 
filter are given without proof. References £lj and [3^ describe 
proofs and [4] gives a detailed discussion for a general case. 

The optimum estimate of the plant states is given by the smoothing 
equation: 

x(k/k) = x(k/k-l) + G(k) [(z(k) - z(k/k-l)J (6-1) 
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0 



where 

x(k/k-l) - $x(k-l/k-l) + AE Wpi'A) 
is the predicted value of the states based upon the last best 
estimate of the states. The expected value of the environmental 
noise is zero because it has been assumed to have a normal 
distribution with zero mean. 

z.(k) = y_(k) + V(k) is the current value of the noisy 
observable. 

a r 

is the predicted noisy 



z(k/k-l) = H x(k/k-l) + E 






observable based on the last value. The expected value of 
measurement noise is zero because it has been assumed to have 
a normal distribution with zero mean, 
and 



^(k) = H x(k) is the current non-noisy value of the obser- 
vable states . 

The filter design problem is based upon the proper selection 
of the gain matrix, G(k), so that the sum of the variances or 
errors associated with the estimate of individual states is 
minimized. The following two recursive relationships are re- 
quired to calculate the filter gain matrix and the conditional 
covariance matrix: 

G(k) - |~P (k/k-l)H T (HP(k/k-l)H T + R) J _1 (6-2) 

and 

P(k/k-l) = <f> [p(k-l/k-2) - G (k- 1) HP (k- 1/k- 2) J $ T + Q 

(6-3) 

where 

H = (1 0) is the observability matrix 
R = E £v*V t J is the covariance matrix of measurement 
noise. Since there is only one observable in the present problem 
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R is a scalar. 

Q = A E [w*W T J A'*' is the covariance matrix of excitation 
or signal noise. A discrete schematic of the filter is shown in 
Figure 1-1 . 

A stable gain matrix of 



G = 



0.786 

0.682 



was obtained from the Fortran 60 program 'Filter' by providing 
it with the following parameters: 



2 

O’ = 0.052 
w 

R(1 , 1) = 0.02 



and 



P(O) = 



0.02 

0 



0 

1 



A listing of program ’Filter' is given in Appendix II. Table 
II— 2 defines the symbolic terms for this program. 

7 . The Filter-Controller 

The combined filter-controller simulation is shown in 
Figure 1-1 . The double lines represent the transmission of 
vectors and the single lines, scalars. Deterministic inputs 
are introduced and are represented by a vector quantity, Dl(k) , 
with the elements consisting of the input and its successive 
derivatives. The difference vector, _B(k), between the determ- 
inistic inputs and the predicted states is multiplied by the 
optimum controller gain matrix, A/*- . A measure of the perform - 
ance of the system can be deduced from the B vector. When 
the B vector is very near zero, the filter is estimating states 
which are approximately equal to the deterministic inputs. 



Then the plant control becomes very small and the system has 
responded to the inputs. 

The output of the controller is a scalar control u(k) which 
is constant over the sampling interval and is given by 

u(k) = AT £x(k) - DI_(k)J (7-1) 

This control must be introduced so that it identically affects 
both the plant and the filter. This is done by feeding the same 
control to the filter as the plant but delayed by one sampling 
interval to account for the difference in the sample indicies. 

The controlled plant and the controlled filter equations are 
then given by 

x(k+l) = $ x(k) + Au(k) + Aw(k) (7-2) 

and 

x(k) = $ x(^"l) + A u(k-l) (7-3) 

Since the train plant of the gun is a type one system, it 
has a steady state error in response to a ramp input. This means 
that the ji vector would not go to zero unless the DI_ vector is 
modified. This is accomplished by introducing a DA factor. 

In this case study, the DI vector is made up of two elements: 
DI(1,1) = step + ramp*t + stabilization 
DI(2 ,1) = ramp 
where 

step is the magnitude of the step 
ramp is the magnitude of the ramp 

and 

stabilization is the magnitude of the stabilization signal 
over one sample interval . 

Consider Figure 7-1 and a unit ramp input. 
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DI 1 = t 



Xj (00) = t 

x 2 (00) = 1 di 2 = 1 

x 2 ( 00 ) = 0 

U(00) = 0 

but 

X 2 (00) = 3U(00) - 3.5X 2 (00) 

0 / -3.5 

hence introduce DA factor 

X 2 (00) = 3 DA - 3.5X 2 (00) 

DA = 3 . 5/3 

A general expression for the DA factor for any magnitude ramp is 
given by 

DA = 3 . 5 ramp ^ ^ 

3 

A Fortran 60 program, 'Hybrid II 1 , was used to simulate the 
filter-controller design and to compare the results with the hy- 
brid simulation which is described in the following sections. A 
listing of this program is included in Appendix II. 

8 . • Description of Hybrid Simulation Equipment 

The experimental part of this thesis was carried out in the 
Hybrid Control Laboratory; United States Naval Postgraduate 
School. Simulation of the system required the following equip- 
ment: 

A. The Control Data Corporation 160 digital computer is 
an electronic computer controlled by an internally stored program 
in sequential locations. Memory capacity is 4096 12-bit words 
of magnetic core storage, with a storage cycle time of 6.4 
microseconds. Instructions are executed in one to four storage 
or memory cycles with the time required for execution varying 
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from 6.4 to 25.6 microseconds or approximately 15 microseconds 
on the average. A computer word is made up of 12 binary digits 
numbered from right to left from 0 to 11 . All positive numbers 
must have a zero in bit 1 1 ; all negative numbers , a one in bit 
11 . Hence positive numbers range from (0000)g to (3777)g and 
negative numbers from (7777)g to (40 0 0) g . 

The CDC 160 is very limited in its arithmetic computation 
capability. Twelve bit addition and subtraction is accomplished 
in two or three memory cycles using ones complement notation. 
Multiplication and division can only be accomplished by success- 
ive 12 bit addition or subtraction, respectively, resulting in 
excessive programming and execution time requirements. 

B. The Control Data Corporation 168 arithmetic unit pro- 
vides the CDC 160 computer with single or double precision, 
fixed point arithmetic operations, either integer or fraction but 
not floating point, for addition, subtraction, multiplication and 
division. However, library subroutine mulop provides single 
precision multiplication in floating point format by first mul- 
tiplying the numbers and then dividing by a scale factor which 
in effect post-shifts the answer so that the binary point is in 
the correct position. Unfortunately, the execution time for 
mulop is very long; 400 microseconds for no scaling and 550 
microseconds for scaling. 

C. The Pace TR20 analog computer is a solid state computer 
having 20 operational amplifiers with a saturation voltage of 
plus or minus ten volts. By the connection of the operate-reset 
relays of the TR20 to the same relays of the ADC-DAC, remote 
control of the analog computer is possible. The run key of the 
CDC 160 starts the digital program and operates the analog 
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computer; the master clear key resets the analog computer. 

D. The ADC-DAC is a 12 bit conversion unit that functions 
in the voltage range of zero to minus ten and has a common minus 
five volts bias source. There are 12 channels analog-to-digital 
and four channels digital-to-analog. Approximately 100 micro- 
seconds is required for each A/D number conversion and approx- 
imately 20 microseconds for each D/A number conversion. The 
input disable jack allows external timing control of the ADC 
channels . 

E. The EAI 1110 Varipl otter is a solid state x-y pen 
recorder with a sensitivity of 100 microvolts per inch. It was 
used to record the simulation results. 

9 . Signal Conditioning and Number Scaling 

Since the voltage ranges of the analog computer and the 
ADC.-DAC, and the number range of the CDC 160 computer are 
not compatible, a certain amount of signal conditioning and 
number scaling is required. Table 9-1 shows the ranges for 
each piece of equipment and their equivalent values . An 
Ax + B transformation on the analog voltages is performed to get 
them into the range of the ADC and thence into the digital com- 
puter. A reciprocal transformation is performed for the transfer 
in the other direction. The A part of the transformation is a 
scale factor equal to 5/8 for the output analog voltages and 
equal to 8/5 for the input analog voltages. This division of 
the output analog voltages by 8 or multiplication by 2"^ pro- 
vides a shifting of the binary point 3 places to the right. Hence 
all voltages being transferred to the CDC 1 60 computer have a 
scale factor of three binary. The B part of the transformation is 
a five volt bias factor required to convert all analog voltages 
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COMPARISON OF EQUIPMENT RANGES 






AND EQUIVALENT VALUES 




ANALOG 


SCALED ADC OUTPUT 


CORRECT 


OBSERVED 


VOLTAGE 


ANALOG - 


- (0 . 625X+Bias) 


DI GITAL 


DIGITAL 


(DECIMAL). 


VOLTAGE 




VALUE 


VALUE 


X 


. 625X 




(OCTAL) 


(OCTAL) 


-8 


-5.000 


- 0.000 


4000 


4003 


-7 


-4.375 


-0.625 


4'377 


4407 


-6 


-3.75 


-1.250 


4777 


5007 


-5 


-3.125 


-1.875 


5377 


5363 


-4 


-2.500 


-2.500 


5777 


5773 


-3 


-1.875 


-3.125 


6377 


6363 


-2 


-1.250 


-3.750 


6777 


7003 


-1 


-0.625 


-4.375 


7377 


7403 


0 


0.000 


-5.000 


0000/7777 


0003 


1 


0.625 


-5.625 


0400 


0377 


2 


1.25 


-6.250 


1000 


0777 


3 


1.875 


-6.875 


1400 


1377 


4 


2.500 


-7.500 


2000 


1777 


5 


3.125 


-8.125 


2400 


2377 


6 


3.750 


-8.750 


3000 


2777 


7 


4.375 


-9.375 


3400 


3377 


8 


5.000 


10.000 


3777 


3773 


TABLE 9-1 




NUMBER SCALING 






DECIMAL 


OCTAL 


BINARY 


SCALED OCTAL 








COMPUTER WORD 


1.000 


1.000 


s 00/1 00/0 00/0 00/0 


0400 


-1.203 


-1.151 


sOO/1 00/1 10/1 00/1 


7313 


-.343 


-.257 


sOO/OOl/Ol 0/1 11/1 


7650 


.020 


.012 


sOO/OOO/OOO/lOl/O 


0005 


.004 


.002 


sOO/OOO/OOO/OOl/O 


0001 



TABLE 9-2 
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into the range of the ADC. 



By limiting the analog voltages on the TR20 to plus or minus 
eight, taking 0.625 of each variable to be outputed to the ADC, 
and adding it to a bias of plus five volts, the output of the 
summer is in the range of zero to minus ten volts ready for in- 
put to the ADC. To go in the other direction, five volts are 
added to the DAC output and the entire quantity is multiplied <by 
1.6. Then the output of the amplifier is in the same range of the 
analog computer. See Figure 11-1 showing this signal transfor- 
mation for program Optimum Controller. 

Since bit number 11 (the 12th bit of the computer word) is 
the sign bit, then the remaining 11 bits are available to repre- 
sent a number. If the decimal point of the computer word is 
scaled three places to the right, this gives the following bit 
representation: 



This implies that numbers in the range of ±7.999 with a binary 
scale factor of three can be accomodated by the computer. 

Since a scale factor of three was selected for numbers in 
the digital computer, all constants must be changed to the 
above format. Table 9-2 shows the decimal, octal, binary and 
scaled computer word representations for a few constants . Sub- 
routine mulop has a second argument called the scale factor. 

In multiplying two numbers with scale factors of three, the 
result has a scale factor of six. Hence a three bit left shift 
of the binary point must be performed by dividing by the scale 



s: 

b] 




xxx 



xxx 
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factor. 

Since the sign bit must be included one bit of the number is 
lost and, therefore, the smallest decimal number that can be 
represented with a binary scale factor of three is .004. However, 
because of inaccuracies in the ADC-DAC conversion as is shown 
in Table 9-1, the smallest decimal number used in the hybrid 
simulation was .02. 

10. Arithmetic & Matrix Subroutines 

Subroutine mulop was the only available arithmetic sub- 
routine at the start of this study. Therefore, before any sim- 
ulation could be attempted, basic arithmetic and matrix sub- 
routines had to be written. These subroutines are described in 
detail in Appendix VI. Error halts are provided in these sub- 
routines for both summation and product overflow. However, 
no provision is made for over-riding this overflow and using the 
full positive or negative values . These subroutines are written 
for single precision or 12 bit numbers. The scaling of the num- 
bers in the computer is flexible and can be changed by modifying 
all constants entered into the computer. The scale factor in 
cell (0027)g must be adjusted, as mentioned in Section 9, to 
return the binary point to the same place after multiplication or 
division . 

1 1 . Program Optimum Controller 

To get a feel for the hybrid simulation problem and to check 
out the arithmetic and matrix subroutines it was decided to try 
a simulation that was simple and whose results could be easily 
verified. Such a simulation problem is the optimal discrete- 
time control system described by the following equation: 
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u(k) =A T [x(k) - DI(k)] (10-1) 

where the true state vector X(k) is assumed completely observ- 
able without noise contamination. This control law, using feed- 
back coefficients obtained from program’ 'Optcon' for the minimum 
time case, will drive the state variables to zero for the regulator 
problem of initial conditions and to some specific value for the 
controller problem of deterministic inputs. The system sim- 
ulation schematic is shown in Figure 11-1 . The listing and 
method of operation of program Optimum Controller are described 
in Appendix III. 

12 . Simulation Results for Program Optimum Controller 
Figures 12-1 through 12-7 are the x-y recorder plots of 

the gun train plant response for various initial conditions and 
deterministic inputs that were limited to the optimal control 
region of the phase plane. In all cases optimum control is 
* effectively achieved in the minimum time of two samples . The 
finite response time indicated on the control signal u(t) is due 
to the inertial lag of the EAI Variplotter. The external timing 
control was found to be the most accurate, and was used for all 
recordings. There were no observable hybrid control problems 
with this relatively simple second order plant and limited inputs. 
However, for higher order plants, hybrid limitations may be 
encountered as discussed in the next section. 

The approximation of the roll and pitch stabilization sig- 
nals by the output of a zero order hold proved to be satisfactory. 
This approximation is used in the filter-controller simulation of 
Section 1 8 . 

13. Hybrid Control System Limitations 

The response of a system controlled by a digital computer 
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FIGURE 12-2 
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FIGURE 12-3 
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FIGURE 12»4 
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FIGURE 13-7 
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may be affected by the following general limitations of hybrid- 
ization: 

A . There is a time delay due to the serial conversion of 
analog quantities to digital numbers by the ADC. The serious- 
ness of this problem is determined by the number of simultaneous 
samples required. In program Optimum Controller, where four 
A/D conversions per sampling interval are required for this 
second order plant, this difficulty was not experienced. If, 
however, the order of the plant was increased, then this effect 
might be noticeable. In such a case an external sample and 
hold device would be necessary to sample all inputs at the same 
instant of time and hold these values until the ADC completes 
the conversion on all channels . 

B. There is a time delay due to the serial nature of the 
computation of the control in a digital computer. Because of 
this step-by-step performance of arithmetic operations, a 
finite time exists between input to the digital computer of 
samples from the plant and output of the control from the com- 
puter to the plant. This time delay is a function of computer 
speed and the complexity of the computation. This problem can 
be minimized by using a high speed digital computer with faster 
arithmetic subroutines than the CDC 160 and by efficient com- 
puter programming between the sampling and the output. 

C. There is a magnitude limitation imposed hy the A/D 
or D/A conversion processes. This is a physical equipment 
limitation in this case with a range of plus or minus five volts 
as discussed in Section 9. In general this magnitude limitation 
will depend upon the specific ADC-DAC equipment, This prob- 
lem would be greatly reduced with shaft to digital encoders , 
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for example. 

D. There is a limit in the magnitude of the floating point 
constants imposed by the computer word length. As discussed in 
Section 9, the smallest number that can be represented in the 
computer with a scale factor of three is .004. A computer with 
an 1 8 bit word length would allow more flexibility in the com- 
putation and entering of constants . 

E. There is a sample timing control problem in that the 
computer is required to sample the inputs at regular intervals . 
This can be done either internally by subroutine delay or ex- 
ternally by a clock. The internal time delay is not as accurate 
as the clock because the actual time delay depends upon vari- 
ations in computation time of arithmetic operations that may 
vary from cycle to cycle. This fact also makes the internal time 
delay much harder to adjust. External timing control was used 
in the hybrid simulation schemes. This type of external clock 
control is much more accurate and would allow time sharing of 
the computer with other equipments . 

14. Selection of Sampling Interval 

In the selection of a sampling interval for the digital 
filter-controller a number of factors must be considered . The 
inertia of the gun is large and the system is velocity and 
acceleration limited. Too small a sampling interval would re- 
quire controls of such large magnitude in the minimum time 
case that they may cause the system to saturate and operate 
non-linearly . On the other hand, for too large a sampling 
interval the late updating of the plant may be too slow for a 
fast target. If a number of systems are to be time shared by 
the same computer, then a long sampling interval would allow 
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more systems to be controlled by the same computer. The types 
of input signals must also be considered. For instance, in 
this case the sampling interval must be small enough in order 
that the roll and pitch stabilization signals can be well rep- 
resented by the output of a zero order hold. The limitations of 
the equipment to be used must also be considered. In this case, 
the parameters and constants had to be within the number range 
of the computer. Considering all these factors, a sampling 
interval of one second was chosen. 

15 . Selection of Noise Values 

The choice of the values of measurement arri environmental 
noise, and the initialization of the error covariance matrix can 
be by conjecture, actual measurement of the statistics or by 
equipment limitations. In this simulation, in order to stay with- 
in the number range of the computer, the selection was predi- 
cated by the latter method. It is assumed that the probability 
density functions of V(k) and W(k) are normal with zero mean. 
Both V(k) and W(k) are one dimensional random variables. Hence 
they are scalars. The noise values for the hybrid simulation 
were selected as follows: 

A. The measurement noise is given by 

R (1 , 1) = E [v(k) 2 ] = 0.02 

The standard deviation of the measurement noise is given by 
o v = 0.1414 

If a shaft to digital encoder is used to transform the shaft angle 
output to a digital number, then the ambiguity of this transfor- 
mation determines the value of R(1 ,1) . 

B. The covariance matrix of random disturbance due to 



44 



environmental or excitation noise is given by 
Q = A E ^W(k) j i = A a ^ A* 



where 



and 



o = W(k), is the variance of excitation noise, 
w 



A A 



t _ 



0.384 0.515 

0.515 0.690 

The variance of excitation noise is calculated by assuming a 
pseudo noise to signal ratio of R(1 ,1)/Q(1 ,1) = 1. 

Therefore 

Q (1 , 1) = .02 = a 2 A A t 
w 

From this relationship the variance and the standard deviation 
of the excitation noise are calculated to be 

a 2 = .052 
w 



and 



<r = .228 
w 

For this weapon, the excitation noise can be considered to be 
a random signal due to jitter. The jitter may be caused by 
noisy transmission of gun orders, sharp variations in the sta- 
bilization or designation signals , and ambiguity in bearing 
alignment between the designator and the gun. 

C. The estimation covariance matrix for the second order 
plant is a two-by-two matrix. P q (1 ,1) is the initial variance 
on position and is based on the statistics of the measurement 
noise. It is therefore reasonable to assign R(1 ,1) as the value 
of P Q (1,1). P Q (2,2) is a measure of the confidence in the 
initial velocity measurement. In this simulation, the velocity 
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of the plant is unknown and must be estimated by the filter. 
Consequently, a large value compared to P Q (1 ,1) is chosen to 
show this uncertainty. Off-dimensional terms are assigned zero 
value since the initial state estimates are uncorrelated. The 
following error covariance matrix was used in program Hybrid II 
to calculate filter gains on-line: 



P 



o 



0.02 

0 



0 

1 



Figure 15-1 represents a normalized plot of the gain matrix 
elements for various pseudo noise-to- signal ratios. R(1 ,1) is 
a measure of the expected position measurement noise power 
and Q(l,l) is a measure of the expected position signal power. 
This plot is made by considering a pseudo noise-to- signal ratio 
in which R(1 ,1) is made equal to multiple values of Q(1 ,1). For 
a noise to signal ratio of one, the stable gain matrix elements 
are always given by 

G(1 ,1) = 0.68 



and 



G (2 , 1 ) = 0.45 
16. Programs Hybrid I & II 

Programs Hybrid I and II were written in machine language 
for the CDC 160 computer to simulate the discrete digital filter- 
controller part of Figure 1-1 . Hybrid I is described thoroughly 
in Appendix IV, and uses constant filter and controller gains. 
Hybrid II is described in Appendix V, and uses on-line com- 
putation of filter gains. These programs solve the controlled 
plant equation 7-2 and the controlled filter equation 7-3 to 
calculate the digital control for the weapon. Figure 16-1 shows 
the simulation schematic for the hybrid control system. 
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In the writing of the Hybrid programs every effort was made 
to minimize the hybridization limitations discussed in Section 13. 
For example, in the first edition of Hybrid II, the time between 
sample and output was 250 miUiseconds, However, by the 
application of good programming techniques and by the re- 
location of certain computation segments, the delay was reduced 
to 130 milliseconds. 

Figure 16-2 shows the time history for programs Optimum 
Controller, Hybrid I and II. In all three programs the A/D 
sampling takes 820 microseconds or less and therefore this 
limitation is small compared to the computation delay times . A 
computation delay of 13.6 milliseconds did not affect the res- 
ponse of the system as observed and discussed in Section 12 for 
the simulation results of program Optimum Controller. However, 
a computation delay time of 130 to 150 milliseconds for the 
Hybrid programs does affect the response of the system as shown 
in Figures 16-3, 18-1 and 18-2. In these Figures a comparison 
of the step response of the system is made between normal gun 
operation being sampled every second, and the system time 
scaled by a factor of ten being sampled every ten seconds. 

Figure 16-4 illustrates in more detail what happens when this 

compute time 't 1 is significant. The time axis of this Figure 
c 

is expanded. If t_ = 0, then the first and succeeding controls 
are based on the samples of the plant at the proper instant and 
the controls are applied to the plant at that same instant pro- 
ducing the desired optimum response. However, if t >0, then 

c 

the control is applied t c seconds after each sample time. In 
the mean time the plant states are advancing without the proper 
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control. When the control is finally applied it is no longer the 



optimal control for the plant at time KT + t . 




O 

N 

EFFECT OF COMPUTATION TIME DELAY 
FIGURE 16-4 

The computation times for the hybrid programs are about 
ten times longer than that for program Optimum Controller. The 
Hybrid programs have an increased complexity of computation but 
the predominant reason for this longer time delay is that a very 
slow multiply subroutine is called very frequently in the matrix 
calculations. Mulop takes 550 microseconds each time it is 
called. It is felt that with newer, higher speed computers' 
and faster arithmetic subroutines, an improvement in computation 
time of ten to one is feasible. Therefore the results obtained by 
time scaling the system by a factor of ten are justifiable rep- 
resentation of the type of control that could be achieved with 
a high speed computer. Because of this consideration, the 
remainder of the simulations results for the Hybrid programs 
are for the time- scaled plant. 
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17 . Adaptive Filter Gains 

In program Hybrid I the steady state filter gains are entered 
as constants. In program Hybrid II the filter gains are computed 
on line for each sample by calling subroutine vf gain. In this 
adaptive filter gain case, the filter gains increase from initially 
computed values, to steady state values in approximately 
seven samples. From this instant on, the adaptive filter gains 
have reached their steady state values and both Hybrid programs 
have the same response for this time -invariant plant. The sim- 
ulation results for programs Hybrid I and II indicate that there 
is no improvement in the response by using adaptive filter gains . 
In fact duplicate pen recordings of the phase plane were taken 
for step responses with one response traced on top of the other. 
Because of this exact reproduction of responses all the simul- 
ation results shown are from one program. Hybrid II. 

1 8 . Simulation Results for Hybrid II 

Figures 18-1 through 18-6 are pen recorder graphs of the 
gun train plant responses for various deterministic inputs of 
step, ramp and stabilization signals, and for measurement and 
environmental noise . Two points must be considered in observ- 
ing the simulation results that contain either noise and/or 
stabilization signal. The noise is reinitialized to start at the 
top of the list for each pen recording so that each individual 
recording can be compared for that particular deterministic 
input. The sinusoidal stabilization signal is not synchronized 
to start at the same position for each pen recording, and there- 
fore, the individual time axes are not correlated. 

Figures 18-1 and 18-2 illustrate again the hybridization 
problem that occurs when too long a delay between sample and 
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output of the control exists. In contrast, the time scaled step 
responses are the desired dead-beat, optimum responses. 

Figure 18-3 shows how the addition of environmental and meas- 
urement noise affect the gun response. In Figure 18-4 a pitch 
stabilization signal is introduced with a step input. The approx- 
imation of the pitch stabilization signal by the output of a zero 
order hold is seen to be satisfactory. In Figure 18-5 the 
requirement for the DA factor to ensure zero steady state error 
of this type one plant in response to a ramp is illustrated. 

Figure 18-6 depicts a typical situation of the weapon slewing in 
response to a step to a particular bearing, and tracking a target 
in response to a ramp. The hump is produced in this position 
response because the DA factor is introduced three samples 
before the ramp is :ihitiated. In hind sight, a more judicious 
method of introducing this DA factor would be to detect the 
magnitude of the deterministic velocity input and calculate 
equation 7-4 each sample. Then, if the deterministic velocity 
vector is zero, the DA factor would also be zero. 

1 9 . Advantages of Digital Control 

The use of a digital filter- controller has many significant 
advantages over the conventional analog system. In cases 
where more accuracy and higher sensitivity are required, 
digital compensation techniques excel over analog methods. 

In addition, virtually noise-free transmission of data in the 
form of digital numerical codes is possible. Non-linear pro- 
gramming and self-optimizing programming techniques can be 
utilized with the digital computer to provide a degree of flex- 
ibility not heretofor achievable with continuous control systems . 

High operating speeds and efficient input/output procedures 
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have made the utilization of general purpose digital computers 
in real-time sampled data systems feasible. Time sharing of 
many equipments with the same computer is another important 
aspect. The use of vector-matrix notation for systems that 
can be represented by linear differential equations is ideally 
suited for digital computation, and has led to development of 
programs for analysis and control of systems. 

20. Conclusions 

The following statements are cited as being important con- 
clusions to this hybrid simulation control problem: 

A. The results of this hybrid simulation have demonstrated 
the feasibility and simplification in design which accompany 
the use of a general purpose computer as an integral component 
of the weapon control system. 

B. The results also demonstrate that the application of a 
discrete filter- controller to provide control of a tactical weapon 
is feasible. 

C. The results demonstrated that there was no difference 
between using constant or adaptive filter gains for this time- 
invariant weapon plant. This means that the simpler com- 
putational program of Hybrid I can be used, and thereby 
minimize the computation time required to effectively control 
the gun. 

D. The application of good programming techniques to 
minimize hybridization limitations is important in providing 
optimum digital control. 

E. The time-shared control of a number of weapons sys- 
tems serially with the performance of other functions is feasible 
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for a general purpose, high speed computer having fast input/ 
output capabilities and fast arithmetic subroutines . 
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APPENDIX I 

STATE SPACE REPRESENTATION OF A SYSTEM 



1-1 General Theory 

The state of a system is defined as a set of variables Which 
in conjunction with the system inputs completely describe the 
behaviour of the system. Consider a linear, time invariant 
system described by the second order differential equation: 
y + ay + c = f(t) (I— 1 ) 

The state variables^ and X 2 and the plant control are defined 
as follows: 
xi = y 

x 2 = x i = y d-2) 

u = f(t) - c 

If f(t) is known for all t greater than or equal to t Q and if the 
initial states are known, then the states can be determined for 
all time t. The second order linear differential equation is 
reduced to a set of first order state equations as follows: 

*1 = x 2 

X£ = -b*Xj - a*X 2 + u 
or written in matrix notation as: 

x = Fx + Du Xl-3) 

where 

x is an nxl vector of the system states 
u is an rxl vector of the system inputs or controls 
F is an nxn matrix of constants 
D is an nxr matrix of constants 
here and 



To f 


D = 


0 


|-b -a 




1 
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In general the output of the system may be a linear combin- 
ation of the system states and controls, and a new set of output 
equations can be defined: 

% = A x + B _u 
where 

y is a qxl vector of system outputs 
A is a qxn matrix of constants 
B is a qxr matrix of constants 

It can be verified that the continuous solution to equation 
Irl is given by 

x(t) = $(t-t 0 ) x(t) + 5 ^(t-t x ) D u(ti) dt 1 (1-4) 

where 

$(t-t Q ) is called the state transition matrix and is defined 
as the inverse Laplace transform of (sI-F) - * 
and 

5 #(t-t 0 ) D dt-^ is called the state distribution matrix 

to 

for the case when the control u is the output of a zero order hold. 
The discrete solution to equation 1-4 is given by 
x(k+l) = $ x(k) + A u(k) (1-5) 

where 

*= c FT 

A = j c F ! t ' tl ! D dt x 
o 

T = t-t Q , T is the sampling interval 

and 

t Q =kT, t = (k+l)T 

The advantages of describing an nth order linear differential 
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equation by a set of first order differential equations are: 

(i) ease of solution of these first order equations by digital 
computer programming techniques 

(ii) ease of correlation between state variables, flow graphs 
and analog computer simulation allowing unified study of sampled- 
data systems 

'(in) ease of specification of initial conditions. 

1-2 . Flow Graphs 

The flow graph of a linear, time invariant system can be 
obtained from the state equation 1-3 or vice versa. Figure 1-1 
was drawn with the following ground rules: 

(i) The system states are selected as outputs of integrators 
and the nodes selected as states have only one branch entering 
the node. 

(ii) x. is connected to x._^ by a unity transmission branch 

(iii) all feedback paths are from the integrator outputs of 
x l/ x 2 etc; to the node representing xj . 
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APPENDIX II 
FORTRAN 60 PROGRAMS 
II- 1 . General Description 

This Appendix contains the following Fortran 60 programs 
that were used to provide constants required for the hybrid 
simulation and also provide a means of verifying the results from 
the hybrid simulation: 

(i) Optcon 

(ii) Filter 

(iii) Hybrid II 

Each program has comment cards which explain their oper- 
ation and a sample data deck is included. Table II- 1 defines 
the terms for program 'Optcon' and outlines the three cases for 
selecting a desired performance index. Table II- 2 defines the 
terms for program 'Filter'. 
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DEFINITION OF TERMS FOR PROGRAM OPTCON 



Variable 


Definition 


Remarks 


N stage 
N 


number of stages 
order of system 




R 


control effort weighting 
constant in the cost 
function 


selected by the operator 
for desired performance 


Q 


state weighting matrix 
in the cost function 


selected by the operator 
for desired performance 


A 


nxl control distribution 
vector 


constant for a given time- 
invariant plant 


$ 


nxn state transition 
matrix 


constant for a given time- 
invariant plant 




intermediate variable 
matrix in recursive 
solution for controller 
gains 


initial value = 0 

o 


P 


intermediate variable 
matrix in recursive 
solutions for controller 
gains 


initial value P Q = Q 


a t 


feedback controller 
gain vector 


initial value A^ = 0 


Case I 


Minimum time or mini- 
mization of terminal states 


Q 0 = I, Q* = °/ R = ° in 

program 'Optcon' 


For second order plant the minimum time performance index is 

J(N) = x 2 (N) + x 2 (N) 

1 ^ 


Case II 


Minimum effort or fuel 


Q o =1, 0= 0, R = 1 in 
program 'Optcon' 


For second order plant the minimum control effort performance 
index is 


J(N) 


= x 2 (n) + x 2 (N) + min 

2 u(k) k=l 


u 2 (k-l) 



TABLE II- 1 
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Case III 



Minimum time and effort Q q = I, Q* ^ 0, R = 1 

in program 'Optcon' 

For second order plant the minimum time and control effort 
performance index is 

J(N) = min ^ x T fc)x(k) + u^(k-l) 

u(k) k=l 

Note: Q and R may be selected to satisfy different cost func- 
tions. Another method for selecting Qand R is given in reference 

[4] • 



TABLE II- 1 (cont'd) 
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12, 12), PS 1(12, 12). P(12, 12) .DEL ( 12 ) ,AT( 12) , 

1 GM CT2 , 12~)TQTY2Tr21TXT 900 ) , I T I T L E ( 1 2) , FM ( 1 ‘ 6TTEM Q21TH MrT2Tf 2T, 

1 Yl (900 ) .Y2 (900 ) .Y3( 900 ) 

C mTS^P^OGRW'CTTTirrZFS^—COSr TUNCT I ONT~J ( N ) =M I NT MUM (‘SUM - XTNTT *Q *XTN')T 

C SUM R*U(N-1)**2) . AN UNLIMITED NUMBER OF ITERATIONS MAY BE MADE AT 
C A COMPUTATION RATE OF 2000 PER MINUTE AFTER THE PROGRAM HAS BEEN 

C COMPUTED". IXE - OUTPUT OF TH I S~ PROGRAM' - ! "S~T H E~ F E E D3"ACK~GATN~TMAT R TX» 

, C A TRANSPOSE. THE FOLLOWING RECURSIVE EQUATIONS WERE DERIVED USING 

C DYNAMIC PROGRAMMING, 

7 — c ATTTCr=~nOErr* PIK-TT^PHTT/TDEITT* PTK=TT* D EU-FRT IT) 

C PSI <K) =PHI+DEL*AT(K> (2) 

C P(K)=PSIT(K)*P(K-1)*PSI (<)+Q+R*A(K)*AT(K) (3) 

— c ptot= otattoi =cr,— psrror=o 

C EQUATIONS 1, 2. AND 3 CONSTITUTE THIS PROGRAM 

C CALX I NDATA - ATTD I'NTTTAtn ZE 

DO 1111 1=1,3 
READ 30.KASE 

3-C-FORMA-T (IT) = 

READ 1 .N.NSTAGE 

1 FORMAT (8110) 

R E AD-2 ,- R ,-D T 

2 FORMAT ( (4F20.0) ) 

READ 2 ( (Q( I.J) ,J=1,N) ,I=1,N) 

C PR I NT-THE—DATA-RE AD— IN. 

PRINT 100 
100 FORMAT ( 1H1 ) 

PR I NT - " 32 » KASE 

32 FORMAT! 9X»5HKASE= ,13) 

PRINT 3 .N.NSTAGE 

3- F ORMATT//9X ,2 H N= ,T3 ,2 0 X , 7 H N S T A G E = » 1 3 ) 

PRINT 4.R.DT 

4 FORMAT ( /9X,2HR=,F15. 11 ,2X,3HDT=,F 15.11 ) 

PR I NT— 9 . ( ( Q( I-. J ) , J=1 , N ) , I =1 , N ) — 

9 FORMAT ( /9X,7HQ( I , J)=/ (2F15.11 ) ) 

CALL PHIDEL(PHI, DEL, N.DT) 

DO-5— I-l-.N — 

DO 5 J=1.N 
GM( I ,J)=0.0 

' EMTT)'= OtO 

FM ( I ) =0 .0 

5 P{ I ,J)=Q( I , J) 

C TNTFT A trIZ E - P ( 0 ) ~ FO R ~ C A S E~ ' T WO 

I F ( KASE-2 ) 35,36,35 
36 P( 1 ♦ 1 ) = 1 . 0 

PT2 »2T=T.'0 

P( 3.3) =1.0 
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35 CONTINUE 
PRINT 19 


T9 f ORM AT7~772~XT 6HTT STT A GTTT4 X ,7 HAT ( IV 1 ) , 3X » 7 H AT fl ,'27 V4X ♦ 6 HP ( 1 , 1 ) , 
14X * 6HP I 1*2) »4X»6HP(2»1) ,4X,6HP(2,2) ) 

C CALCULATE AT(J) 


DO 27“~1<K=T> WSTA'G'E - 

DEN=0.0 

D06 1 = 1, N 


DOb J= 1 ,N 

6 EMI I ) = EM( I )+DEL(J)*P(J,I ) 
D08 1 = 1, N 


D07 J = 1 , N 

7 FM ( I ) =FM( I ) +EM( J) *PH I I J » I ) 

8 DEN=DEN+EM ( I )*DEL( I ) 


DEN = -DE"N-R 

D010I = 1 ,N 

ATI I ) = FM I I ) /DEN 


FKIT) =0,0 
10 EMI I ) =0 .0 

C CALCULATE PS I I K , I , J ) 


D0T3T=1 *N 
DO 1 3 J= 1 » N 

13 PSI (I,J)=PHI(I,J)+DEL( I ) *AT I J ) 


C~ CA L CUUATE P'f K , I , J ) 

DO 15 1=1, N 
DO 15 J = 1 , N 


DO - IS - L = 1 »N 

15 GM I I , J ) =GM I I , J ) +PS I I L, I )*P (L , J ) 
DO 171=1, N 


DO 17 J = 1 , N ' ' ‘ ' ' 

DO 16 L=1,N 

16 HMI I ,J )=HM( I , J ) +GM I I ,L )*PSI I L, J) 


C CASE 1 TERMINAL CONTROL, OMIT Q(I,J> IN EQUATION FOR PI I , J ) 

C CASE2 MINIMUM FUEL OMIT Q I I , J > IN EQUATION FOR P I I , J > 


I F I KASE-2 ) '31 , 3T , 3 3 
31 P< I,J)=HM(I,J)+R*AT( I)*AT< J) 
GO TO 17 


C CASE THREE MINIMIZATION OF TIME AND FUEL 

3 3 PI I ,J)=HM( I ,J)+Q( I ,J)+R*AT I I )*AT( J) 


— — rruMi i vj) =o.o 

DOl 8 1=1, N 
DO 18 J=1»N 


- - 18 GM I'l ♦ J ) =0 . 0 

PRINT 20»KK,AT I 1 ) , AT I 2 ) ,P I 1 , 1 ) ,P( 1 ,2 ) »P I 2 , 1 ) ,P< 2,2 ) i 

20 FORMAT I 15 , I 6F10.6) ) 


'22 CONTINUE 
1111 CONTINUE 
END 
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SUBROUTINE PH I DEL I PH I , DEL » N , DT ) 

C VALID ONLY FOR A CONSTANT F MATRIX 

DIMENSION- T(T2"riZr.PHT C12 ,T2 )VTERM I 12 .12 T.WORMI Y2TT21 
1 * DEL ( 12 ) »DELM( 12,12) ,TELM( 12 , 12) »DELP( 12.12) *D(12) 

1 FORMAT I (8F10.0) ) 

RSADrrrrFriRTrcrric=TTMT, tr=i ,n ) 

READ 1 (D( I ) , 1 = 1 ,N ) 

1003 PRINT 399»DT»((F(IR»IC)»IC=1»N)»IR=1»N) 

3~9 9~F0 R M A rr/'r/3 H DT=-rrF5'7T/ , 7 H FIT VJ )=/ r '.T(2 F8 Y2TT) 

PRINT 3991 CDtl ),I=1,N) 

3991 FORMAT ( / 5HD I I ) =/ I 2F8 .2 ) ) 

N FTNAX=1 

TM = 0.0 

DO 400 I R=1 *N 

do -4 o o-rc-rrN 

TERMC I R • I C ) =0.0 
WORM! IR.IC) =0.0 

TERM Cl R v I RT=TvO 

TELM ( I R » I C ) =TERM ( IR, IC)*DT 
DELP ( I R . I C ) =TELM ( IR, IC) 

DELM ( I R , I C)-=0v0 — : 

DELI IR)=0.0 

400 PHI ( IR, I C ) =TERM I I R , I C ) 

4 T M = lvO+TM 

DO 500 I R = 1 ,N 
DO 500 I C = 1 ,N 

DO - 500 JN=~rYN ~ 

DELMI IR,IC)=DELM( IR, IO-TELMI IR,JN)*F( JN, I C ) *DT/ ( TM+1,0 ) 
500 WORM! I R ♦ I C ) =TERM( I R, JN )*F ( JN, I C )*DT/TM+WORM( IR.IC) 

D0~401 IR = r,N 

DO 401 I C = 1 , N 

TERM! IR, IC) =WORM( IR, IC) 

TELM M R , I C-)-=DELM M R rIXr) 

DELPI IR, IC) =DELP( IR, I C ) +TELM ( IR.IC) 

PHI ( IR, IC) = PHI I IR, IO+TERMI IR.IC) 

DELM MR ,TCT=0v0 

401 WORM! IR, 10=0.0 
ABC = 0 ♦ 0 

DO 2 1 - 1-, N — 

DO 2 J= 1 ,N 
AA=TERM ( I ,J ) 

AB-ABSF ( AA ) 

IFIABC-AB) 3,3,2 
3 ABC=AB 

2-CON T-INUE 

I FI 0.00000000 5-ABC) 5, 5,6 
5 GO TO 4 

: “6 — PRI NT IT, I I PH I I ITUTTTT=T7 NTVT=1TN7 

11 FORMAT! /9X.9HPH I I I , J ) =/ I2F15* 11 ) ) 

DO 600 1=1, N 
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DO 600 <= 1 » N 
DO 600 J = 1 » N 

wo~ DEL _ rrr= DYrrn^PTTnT7T)*DEL p (ejvkt*dtic)' 

PRINT 1 2 » ( DEL ( I ) » I = 1 > N ) 

12 FORMAT ( 9X,7HDEL ( I ) =/ ( 1 FI 5 . 1 1 ) ) 

end 

END 
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DEFINITION OF TERMS FOR PROGRAM FILTER 



Variable Definition 



Remarks 



A 

G 
* H 



R 

Q 

p 

X 

Y 

\ z 

A 

X 

A 

Z 

* 

X 

W(k) 



nxn state transition matrix 



nxl control distribution vector 



nxl optimum filter gain matrix 
plant observability matrix 



measurement noise covariance 
matrix 

random excitation distribution 
covariance matrix 

error covariance matrix 



nxl vector of plant states 

nxl vector of plant non noisy 
observable states 



constant for a given 
time-invariant plant 

constant for a given 
time- invariant plant 



order of matrix depends 
upon number of plant 
observables 

value selected by 
operator 

value selected by 
operator 

initial value selected 
by operator 



nxl vector of noisy plant ob- 
servable states 



nxl vector of the predicted 
value of X(k), based on the 
previous observation X(k-l) 

nxl vector of the predicted 
value of Z(k), based on the 
previous observation Z(k) 

nxl vector of the best estimate 
of X(k) , based on the current 
observation of Z(k) 

random disturbance of environ- assumed normal with 
ment or excitation noise 



TABLE II— 2 
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mean zero 



Variable 


Definition 


Remarks 


V (k) 


random disturbance of 
measurement noise 


assumed normal with 
mean zero 


0 w 


standard deviation of 
excitation noise 




Ov 


standard deviation of 
measurement noise 




TABLE II-2(cont'd) 
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V. 



PROGRAM F I I..TFR 

C D 1 ORDER OF SYSTEM IN 12 FORMAT 

"C D2~SAM PLTNG I NTERVAin N" 8 F 1 0 . 0 FORMAT 

C D3 F MATRIX 3Y ROWS IN 8F10.0 FORMAT 

C D4 D MATRIX 3Y COLUMN IN 8F10.0 FORMAT 

"C DF5 V'ATRTAA r CE~0'F E'X'CT TATI 0 N _ N 0 1 SE - S I GWSQ - TN - 8 FT0T6 - FORMAT 

C D6 R COVARIANCE MATRIX OF MEASUREMENT NOISE IN 8F10.6 FORMAT 

C PHI SYSTEM TRANSITION MATRIX 

'C OEr~DTSTRT3TfTTON~MATRTX 

C 6 OPTIMUM GAIN MATRIX 

C H OBSERVABLE MATRIX 

'C P~B E ST E STTMATE - 0 F"E R R 0 R~C 0 V A RTA N C E~'M ATRTX 

C 0 EXCITATION NOISE COVARIANCE MATRIX 

DIMENSION P(12,12),Q(12»12),H(l2,12),R(12,12),G(12,12),PHIT(12»12) 
r.PHITl 2TT2 TTD-EEnTTTUE LDELTCT2 VT2 1TDELSTT2TT2T7DEC5T ( YZ7Y21~* 

2 X C 12 ) *XHAT ( 12) »YHAT ( 12 ) »PNEW(l2.12) 

READ 33 » N 

3~3 — FO R M AT (T2T “ 

READ 2 * DT 
2 FORMAT(8F10.0) 

PFTNT - 1 003 

1003 FORMAT ( 1H1 ) 

CALL PHIDELtPHI »DEL»N»DT) 

DO~1001~CL~=T~r4 

READ 20 »S I GWSQ 
20 FORMAT ( F10«6 ) 

D O TOITl - LX =TT4 

C INIALIZE MATRICES FOR EACH RUN 

DO 10 1=1 » N 

D ^5 i' 0 _ GJ=TTN 

PNEW ( I ♦ J) =0.0 
P( I » J ) =0.0 

GT ITUT= 0 VO 

0( I , J) =0.0 
DELDELT ( I » J ) =0 .0 

d e Lsrrvj ) = ovo 

10 DELST ( I » J ) =0 .0 
READ 2001 »R ( 1 * 1 ) 

2 00 1 — FORMAT ( 8F10 .6 ) 

DO 30 1=1, N 
30 DELS ( I » 1 ) =DEL ( I ) 

CAUL TRANS ( DELS', N , N , DELST ) . 

CALL PROD f DELS , DELST ,N ,N,N, DELDELT ) 

DO 40 1=1, N 

D0“4'0~ J = 1 , N — 

40 Q( I ,J)=S I GWSQ* DELDELT ( I ♦ J) 

SIGW=SORTF( SIGWSO) 

PRTNT“r004TS' I GW 

1004 FORMAT ( / 10X,5HSIGW= /,E20.8) 
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